R Notebook

利用{rayshader}生成全球任意位置的三维地形图. refer to https://github.com/camartinezbu/30DayChartChallenge2022/blob/f7f7cb1337f2604da2b78d7c920cb318606cf4c6/14-3-dimensional/scripts/1-plot.R

Load libraries

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ──────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✔ ggplot2 3.3.6     ✔ purrr   0.3.4
✔ tibble  3.1.7     ✔ dplyr   1.0.9
✔ tidyr   1.2.0     ✔ stringr 1.4.0
✔ readr   2.1.2     ✔ forcats 0.5.1
── Conflicts ─────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

library(osmdata)

Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright

library(sf)

Linking to GEOS 3.9.1, GDAL 3.3.2, PROJ 7.2.1; sf_use_s2() is TRUE

library(rayshader)
library(raster)

载入需要的程辑包：sp

载入程辑包：‘raster’

The following object is masked from ‘package:dplyr’:

    select

library(elevatr)

Create Bounding Box

定义地理区域范围.

# map center longitude and latitude
center_lon_lat <- c(102.929, 30.150)
half_width <- 0.05

center_lon_lat <- c(116.321,39.947)
half_width <- 0.1

med_bbox <- st_bbox(c(xmin = center_lon_lat[1]-half_width, xmax = center_lon_lat[1]+half_width, 
                      ymin = center_lon_lat[2]-half_width, ymax = center_lon_lat[2]+half_width),
                      crs = 4326)

med_bbox_df <- data.frame(x = c(center_lon_lat[1]-half_width, center_lon_lat[1]+half_width),
                          y = c(center_lon_lat[2]-half_width, center_lon_lat[2]+half_width))


extent_zoomed <- raster::extent(med_bbox)

Get elevation data

利用{elevatr}下载地形数据.

prj_dd <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

elev_med <- get_elev_raster(med_bbox_df, prj=prj_dd, z=12, clip="bbox")

 Accessing raster elevation [-------------------------]   0%
 Accessing raster elevation [=>-----------------------]   8%
 Accessing raster elevation [===>---------------------]  17%
 Accessing raster elevation [=====>-------------------]  25%
 Accessing raster elevation [=======>-----------------]  33%
 Accessing raster elevation [=========>---------------]  42%
 Accessing raster elevation [===========>-------------]  50%
 Accessing raster elevation [==============>----------]  58%
 Accessing raster elevation [================>--------]  67%
 Accessing raster elevation [==================>------]  75%
 Accessing raster elevation [====================>----]  83%
 Accessing raster elevation [======================>--]  92%
 Accessing raster elevation [=========================] 100%

Mosaicing & Projecting
Clipping DEM to bbox
Note: Elevation units are in meters.

elev_med_mat <- raster_to_matrix(elev_med)

[1] "Dimensions of matrix are: 1359x1358"

Get OpenStreetMap(OSM) data

从OpenStreetMap获得地理信息数据

# Get highway data
# https://wiki.openstreetmap.org/wiki/Map_features#Highway
med_highway <- med_bbox %>%  opq() %>%
  add_osm_feature("highway",
    c("motorway", "trunk", "primary", "secondary", "tertiary")) %>%
  osmdata_sf()
med_highway_lines <- med_highway$osm_lines

# Get river data
# https://wiki.openstreetmap.org/wiki/Map_features#Waterway
med_rivers <- med_bbox %>% opq() %>%
  add_osm_feature("waterway",
    c("river", "stream")) %>%
  osmdata_sf()
med_river_lines <- med_rivers$osm_lines

# Get buildings data
med_buildings <- med_bbox %>% opq() %>%
  add_osm_feature("building") %>%
  osmdata_sf()
med_building_polygons <- med_buildings$osm_polygons

# Get parkings data
med_parkdings <- med_bbox %>% opq() %>%
  add_osm_feature("parking") %>%
  osmdata_sf()
med_parking_polygons <- med_parkdings$osm_polygons

# Get tourism data
med_tourisms <- med_bbox %>% opq() %>%
  add_osm_feature("tourism") %>%
  osmdata_sf()
med_tourism_polygons <- med_tourisms$osm_polygons

# Get place data
med_places <- med_bbox %>% opq() %>%
  add_osm_feature("place") %>%
  osmdata_sf()
med_places_points <- med_places$osm_points %>%
  filter(place %in% c("town","city"))

Draw base map

base_map <- elev_med_mat %>% 
  height_shade() %>%
  add_overlay(sphere_shade(elev_med_mat, texture = "desert", colorintensity = 5), alphalayer=0.5) %>%
  add_shadow(lamb_shade(elev_med_mat), 0) %>%
  add_shadow(ambient_shade(elev_med_mat),0) %>% 
  add_shadow(texture_shade(elev_med_mat,detail=8/10,contrast=9,brightness = 11), 0.1)

plot_map(base_map)

Add OpenStreetMap(OSM) data

# create highway layer
med_road_01 <- med_highway_lines %>%
  filter(highway %in% c("motorway","trunk"))
med_road_02 <- med_highway_lines %>%
  filter(highway %in% c("primary", "secondary"))
med_road_03 <- med_highway_lines %>%
  filter(highway %in% c("tertiary"))

road_layer <- generate_line_overlay(med_road_03,extent = extent_zoomed,
                                    linewidth = 2, color="white", 
                                    heightmap = elev_med_mat) %>% 
  add_overlay(generate_line_overlay(med_road_02,extent = extent_zoomed,
                                    linewidth = 1, color="white", lty=3,
                                    heightmap = elev_med_mat)) %>% 
  add_overlay(generate_line_overlay(med_road_01,extent = extent_zoomed,
                                    linewidth = 3, color="white",
                                    heightmap = elev_med_mat))

# create wateray layer
med_river_01 <- med_river_lines %>%
  filter(waterway %in% c("river"))
med_river_02 <- med_river_lines %>%
  filter(waterway %in% c("stream"))

water_layer <- generate_line_overlay(med_river_01,extent = extent_zoomed,
                                     linewidth = 2, color="skyblue2", 
                                     heightmap = elev_med_mat) %>%
  add_overlay(generate_line_overlay(med_river_02,extent = extent_zoomed,
                                    linewidth = 1, color="skyblue2", lty=3,
                                    heightmap = elev_med_mat))

# polygon_layer = generate_polygon_overlay(med_parking_polygons, extent = extent_zoomed,
#                                          heightmap = elev_med_mat, palette="grey30") %>%
#   add_overlay(generate_polygon_overlay(med_building_polygons, extent = extent_zoomed,
#                                        heightmap = elev_med_mat, palette="darkred")) %>% 
#   add_overlay(generate_polygon_overlay(med_tourism_polygons, extent = extent_zoomed,
#                                        heightmap = elev_med_mat, palette="darkgreen"), alphalayer = 0.6)
polygon_layer = generate_polygon_overlay(med_building_polygons, extent = extent_zoomed,
                                       heightmap = elev_med_mat, palette="darkred")

Error in s2_geography_from_wkb(x, oriented = oriented, check = check) : 
  Evaluation error: Found 1 feature with invalid spherical geometry.
[28669] Loop 0 is not valid: Edge 7 crosses edge 10.

base_map %>% 
#  add_overlay(polygon_layer) %>%
  add_overlay(water_layer, alphalayer = 0.8) %>% 
  add_overlay(road_layer) %>%
  add_overlay(generate_label_overlay(med_places_points, extent = extent_zoomed,
                                     text_size = 1, point_size = 1, color = "white",
                                     halo_color = "black", halo_expand = 10, 
                                     halo_blur = 20, halo_alpha = 0.8,
                                     seed=1, heightmap = elev_med_mat, data_label_column = "name.en")) %>%
  plot_map()

base_map %>% 
#  add_overlay(polygon_layer) %>%
  add_overlay(water_layer, alphalayer = 0.8) %>% 
  add_overlay(road_layer) %>%
  add_overlay(generate_label_overlay(med_places_points, extent = extent_zoomed,
                                     text_size = 1, point_size = 1, color = "white",
                                     halo_color = "black", halo_expand = 10, 
                                     halo_blur = 15, halo_alpha = 0.8,
                                     seed=1, heightmap = elev_med_mat, data_label_column = "name.en")) %>%
  plot_3d(elev_med_mat, windowsize=c(1200,1100), zscale=20, theta=20, phi=30, fov=45, zoom=0.6)